FUNCTION ACA_IIR() provides an estimate of a filter's coefficients using a stochastic gradient method.
Based on the paper of C. Cheong Took and D. Mandic, "Adaptive IIR Filtering of Noncircular Complex Signals", IEEE Transactions on Signal Processing, to be published, 2009.
INPUT:
x: filter input [(N+1) x 1]
d: desired response [(N+1) x 1]
M: number of taps
mu: step-size
OUTPUT:
a, b, h, g : estimated taps [M x 1]
e: estimation error for each tap [1 x n]
y: filter output
y(n)= dot(a,y) + dot(b,x) + dot(g,y^*) + dot(h,x^*) ;
Complex Valued Nonlinear Adaptive Filtering toolbox for MATLAB
Supplementary to the book:
"Complex Valued Nonlinear Adaptive Filters: Noncircularity, Widely Linear and Neural Models"
by Danilo P. Mandic and Vanessa Su Lee Goh
(c) Copyright Danilo P. Mandic 2009
http://www.commsp.ee.ic.ac.uk/~mandic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You can obtain a copy of the GNU General Public License from
http://www.gnu.org/copyleft/gpl.html or by writing to
Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
...........................................0001 % FUNCTION ACA_IIR() provides an estimate of a filter's coefficients using a stochastic gradient method. 0002 % 0003 % Based on the paper of C. Cheong Took and D. Mandic, "Adaptive IIR Filtering of Noncircular Complex Signals", IEEE Transactions on Signal Processing, to be published, 2009. 0004 % 0005 % INPUT: 0006 % x: filter input [(N+1) x 1] 0007 % d: desired response [(N+1) x 1] 0008 % M: number of taps 0009 % mu: step-size 0010 % 0011 % OUTPUT: 0012 % a, b, h, g : estimated taps [M x 1] 0013 % e: estimation error for each tap [1 x n] 0014 % y: filter output 0015 % y(n)= dot(a,y) + dot(b,x) + dot(g,y^*) + dot(h,x^*) ; 0016 % 0017 % 0018 % Complex Valued Nonlinear Adaptive Filtering toolbox for MATLAB 0019 % Supplementary to the book: 0020 % 0021 % "Complex Valued Nonlinear Adaptive Filters: Noncircularity, Widely Linear and Neural Models" 0022 % by Danilo P. Mandic and Vanessa Su Lee Goh 0023 % 0024 % (c) Copyright Danilo P. Mandic 2009 0025 % http://www.commsp.ee.ic.ac.uk/~mandic 0026 % 0027 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0028 % This program is free software; you can redistribute it and/or modify 0029 % it under the terms of the GNU General Public License as published by 0030 % the Free Software Foundation; either version 2 of the License, or 0031 % (at your option) any later version. 0032 % 0033 % This program is distributed in the hope that it will be useful, 0034 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0035 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0036 % GNU General Public License for more details. 0037 % 0038 % You can obtain a copy of the GNU General Public License from 0039 % http://www.gnu.org/copyleft/gpl.html or by writing to 0040 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0041 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0042 % ........................................... 0043 function [a,b,h,g,e,y] = ACA_IIR(x,d,M,mu) 0044 0045 L=length(x); % length of data 0046 N=L-1 ; 0047 StepA=1; % one-step ahead prediction 0048 0049 for stepAhead=1:StepA 0050 0051 y = (zeros(1,L) + i*zeros(1,L))'; 0052 e=y; 0053 0054 a = transpose(zeros(1,M-1) + i*zeros(1,M-1)); 0055 b = transpose(zeros(1,M) + i*zeros(1,M)); 0056 h = transpose(zeros(1,M) + i*zeros(1,M)); 0057 g = transpose(zeros(1,M-1) + i*zeros(1,M-1)); 0058 0059 psia= transpose(zeros(1,M-1) + i*zeros(1,M-1)); 0060 psib = transpose(zeros(1,M) + i*zeros(1,M)); 0061 psig= transpose(zeros(1,M-1) + i*zeros(1,M-1)); 0062 psih = transpose(zeros(1,M) + i*zeros(1,M)); 0063 0064 phia= transpose(zeros(1,M-1) + i*zeros(1,M-1)); 0065 phib = transpose(zeros(1,M) + i*zeros(1,M)); 0066 phig= transpose(zeros(1,M-1) + i*zeros(1,M-1)); 0067 phih = transpose(zeros(1,M) + i*zeros(1,M)); 0068 0069 0070 %begin of algorithm 0071 for n = M : N 0072 if (n+1+stepAhead)>N 0073 break 0074 end 0075 0076 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0077 % INPUTS 0078 u = x(n:-1:n-M+1) ; 0079 conju=conj(u); 0080 0081 if n==M 0082 y(1:M) = u(1:M); 0083 end 0084 0085 yy = y(n-1:-1:n-M+1) ; 0086 conjyy = conj(yy); 0087 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0088 0089 0090 y(n+stepAhead)= dot(a,yy) + dot(b,u) +... 0091 dot(g,conjyy) + dot(h,conju) ; 0092 0093 0094 e(n) = d(n+stepAhead) - y(n+stepAhead) ; 0095 conje = conj(e(n)); 0096 0097 %%%%%%%%%%%%Updates of sensitivities or gradients%%%%%%%%%%%%%%%%%%%%%%%%%% 0098 if n==M 0099 phia = conjyy; 0100 phib = conju; 0101 phig = yy; 0102 phih = u; 0103 0104 else 0105 phia = [conjyy(1) + dot(conj(a),phia) + dot(conj(g),psia); phia(1:end-1)] ; 0106 phib = [conju(1) + dot(conj(a),phib(1:end-1)) + dot(conj(g),psib(1:end-1)); phib(1:end-1)]; 0107 phig = [yy(1) + dot(conj(a),phig) + dot(conj(g),psig); phig(1:end-1)]; 0108 phih = [u(1) + dot(conj(a),phih(1:end-1)) + dot(conj(g),psih(1:end-1)); phih(1:end-1)]; 0109 0110 psia= [dot(a,psia) + dot(g,phia); psia(1:end-1)]; 0111 psib= [dot(a,psib(1:end-1)) + dot(g,phib(1:end-1)); psib(1:end-1)]; 0112 psig= [dot(a,psig) + dot(g,phig); psig(1:end-1)]; 0113 psih= [dot(a,psih(1:end-1)) + dot(g,phih(1:end-1)); psih(1:end-1)]; 0114 end 0115 0116 0117 %%%%%%%%%%%%Updates of coeffs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0118 if n< 101 % only 100 first samples to train, then use same coefficients 0119 mua = mu; 0120 mub = 5*mu; 0121 muh = (10^-1)*mu; 0122 mug = mu; 0123 else 0124 mua=0; 0125 mub=0; 0126 muh = 0; 0127 mug = 0; 0128 end 0129 0130 a = a + mua*( e(n)*phia + psia*conje ); 0131 b = b + mub*( e(n)*phib + psib*conje); 0132 h = h + muh*( e(n)*phih + psih*conje ); 0133 g = g + mug*( e(n)*phig + psig*conje); 0134 end 0135 0136 end 0137 0138 0139 %Plotting of graphs 0140 figure,subplot(211), plot(real(y),'k-.') 0141 hold on,plot(real(x),'r') 0142 legend('ACA IIR','Actual'), grid 0143 axis([1 length(x) min(real(x)) max(real(x))]), xlabel('Time (samples)') 0144 subplot(212), plot(imag(y),'k-.') 0145 hold on,plot(imag(x),'r'), grid 0146 axis([1 length(x) min(imag(x)) max(imag(x))]), xlabel('Time (samples)') 0147 0148